;; Non-robust Ramer-Douglas-Peucker algorithm
;; Copyright (C) 2023  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.ramer-douglas-peucker
  (:require
   [clojure.math :refer [sqrt]]))

;; https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points
(defn- compute-distance [[x1 y1] [x2 y2] [x0 y0]]
  (let [numer (abs (- (* (- x2 x1) (- y1 y0))
                      (* (- x1 x0) (- y2 y1))))
        denom (sqrt (+ (* (- x1 x2) (- x1 x2))
                       (* (- y2 y1) (- y2 y1))))]
    (/ numer denom)))

(defn- dmax-reducer
  [a b [dmax index] [i pt]]
  (let [dist (compute-distance a b pt)]
    (if (> dist dmax)
      [dist i]
      [dmax index])))

(defn- find-dmax
  "Find the element of `points` that is furthest from the line defined by the first and last elements of `points`."
  [points]
  (let [a (first points)
        b (last points)]
    (->> points (map-indexed vector) (reduce #(dmax-reducer a b %1 %2) [0 0]))))

(defn classic-rdp
  "Apply the Ramer-Douglas-Peucker algorithm to `points` using distance `epsilon`."
  [points epsilon]
  (if (< (count points) 3)
    points
    (let [[dmax index] (find-dmax points)]
      (if (> dmax epsilon)
        (let [r0 (classic-rdp (take (inc index) points) epsilon)
              r1 (classic-rdp (drop index points) epsilon)]
          (concat (drop-last r0) r1))
        [(first points) (last points)]))))

(defn heap-rdp
  "Apply the Ramer-Douglas-Peucker algorithm to `points` using distance `epsilon`."
  [points epsilon]
  (if (< (count points) 3)
    points
    (loop [work (vector points)
           result []]
      (if-let [pts (peek work)]
        (let [[dmax index] (find-dmax pts)]
          (if (> dmax epsilon)
            (recur (conj (pop work) (drop index pts) (take (inc index) pts)) result)
            (recur (pop work) (conj result (first pts)))))
        (conj result (last points))))))
